home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_10_04 / 1004028a < prev    next >
Encoding:
Text File  |  1992-02-18  |  3.1 KB  |  103 lines
  1. /* Listing 4. */
  2. typedef struct {
  3.   float X1, X2;
  4. } ARG_F_2;            /* vector 2 */
  5. #include "float.h"
  6. #include <errno.h>
  7. #ifndef errno
  8. extern int errno;
  9. #endif
  10. #define MANTBITS (FLT_MANT_DIG -1)
  11. ARG_F_2 powf_2(xi1, y1, xi2, y2)
  12. union fltfmt {
  13.   float flt;
  14.   int iflt;            /* VAX: must change all this */
  15.   struct ffmt {
  16.     unsigned int ex:9;
  17.     unsigned int mant:MANTBITS;
  18.   } fmt;
  19. } xi1, xi2, y1, y2;
  20. {
  21. #define max(i,j) ((i)>(j)?(i):(j))
  22. #define min(i,j) ((i)<(j)?(i):(j))
  23. #if FLT_MANT_DIG != 24
  24. #error "use portable frexp() ldexp() */
  25. #endif
  26. #if  FLT_ROUNDS == 1
  27. #if defined(__STDC__)
  28. /* This works on some non-ANSI compilers */
  29. #define ROUND(x) ((x)>=0?( \
  30.   (x)+1/LDBL_EPSILON)-1/LDBL_EPSILON: \
  31.         ((x)-1/LDBL_EPSILON)+1/LDBL_EPSILON)
  32. #else
  33. #define ROUND(x) ((x)>=0?( \
  34.   (x)+1/LDBL_EPSILON)*1-1/LDBL_EPSILON: \
  35.     ((x)-1/LDBL_EPSILON)*1+1/LDBL_EPSILON)
  36. #endif
  37. #else
  38. #define ROUND(x) ((x)>=0?(int)(x+.5):(int)(x-.5))
  39. #endif
  40.   int mi, mi2, msign;
  41.   double xr, x2, r, r1;
  42.   ARG_F_2 res;
  43. /* Copy 1 */
  44. /* This frexp() operation would be done better after
  45. ** promotion to double
  46. ** but it's not mandatory unless dealing with gradual
  47. ** underflow;
  48. ** it would eliminate most cases of 0 and Inf changing
  49. ** to finite numbers
  50. ** if((xi1.flt=frexp(xi1.flt,&mi))<sqrt(.5)){
  51.   --mi;
  52.   xi1.flt *= 2;
  53.   } */
  54.   mi = ((xi1.iflt & 0x7fffffff) -
  55.     (mi2 = (xi1.fmt.mant < 0x3504f3 ?
  56.             (2 - FLT_MIN_EXP) << MANTBITS :
  57.             (1 - FLT_MIN_EXP) << MANTBITS)))
  58.     >> MANTBITS;
  59.   if (xi1.iflt < 0 | xi2.iflt < 0) errno = EDOM;
  60.   xi1.iflt = mi2 | xi1.fmt.mant;
  61. /* Mult by y distributed to increase parallelism */
  62.   r1 = (xr = (xi1.flt - 1) / (xi1.flt + 1)) * y1.flt;
  63.   x2 = xr * xr;
  64. /* Coefficients determined by Chebyshef fitting
  65. ** double precision is only really needed from here */
  66.   r = y1.flt * (double) mi + r1 * (2.8853904 +
  67.     x2 * (.5958 * x2 + .961588));
  68. /* Msign = (r -= rint(r)) <0 */
  69.   msign = (r -= r1 = ROUND(r)) < 0;
  70.   r *= 125.0718418 + (x2 = r * r);
  71.   x2 = 360.8810526 + 17.3342336 * x2;
  72. /* Xi1.flt = ldexp((x2+r)/(x2-r),(int)r1) */
  73.   xi1.flt = (x2 + r) / (x2 - r);
  74. /* Preferably do this ldexp() operation in double,
  75. ** but it's slower,
  76. ** even though msign can be eliminated;
  77. ** it would always give Inf rather than NaN
  78. ** and would allow use of gradual underflow */
  79.   xi1.iflt += (max(FLT_MIN_EXP - 2 + msign,
  80.         min(FLT_MAX_EXP + msign, (int) r1)) << MANTBITS);
  81.   /* X.fmt.ex+=mi; with limiting to prevent exponent wraparound */
  82.   res.X1 = xi1.flt;
  83. /* Copy 2 */
  84.   mi = ((xi2.iflt & 0x7fffffff) -
  85.     (mi2 = (xi2.fmt.mant < 0x3504f3 ?
  86.             (2 - FLT_MIN_EXP) << MANTBITS :
  87.             (1 - FLT_MIN_EXP) << MANTBITS)))
  88.     >> MANTBITS;
  89.   xi2.iflt = mi2 | xi2.fmt.mant;
  90.   r1 = (xr = (xi2.flt - 1) / (xi2.flt + 1)) * y2.flt;
  91.   r1 *= x2 = xr * xr;
  92.   r = y2.flt * ((double) mi + xr * 2.8853904) +
  93.     r1 * (.5958 * x2 + .961588);
  94.   msign = (r -= r1 = ROUND(r)) < 0;
  95.   r *= 125.0718418 + (x2 = r * r);
  96.   x2 = 360.8810526 + 17.3342336 * x2;
  97.   xi2.flt = (x2 + r) / (x2 - r);
  98.   xi2.iflt += (max(FLT_MIN_EXP - 2 + msign,
  99.         min(FLT_MAX_EXP + msign, (int) r1)) << MANTBITS);
  100.   res.X2 = xi2.flt;
  101.   return res;
  102. }
  103.